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Abstract 

A two-level, global, spectral model using pressure as a vertical 
coordinate is developed. The system of equations describing the model 
is nonlinear and quasi-geostrophic (linear balance) (Lorenz, 1960). A 
moisture budget is calculated in the lower layer only with moist convec- 
tive adjustment between the two layers. The mechanical forcing of 
topography is introduced as a lower boundary vertical velocity. Solar 
forcing is specified assuming a daily mean zenith angle. On land and 
sea ice surfaces a steady state thermal energy equation is solved to 
calculate the surface temperature. Over the oceans the sea surface 
temperatures are prescribed from the climatological average of January. 
The model is integrated to simulate the January climate. 
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1. Introduction 


Numerical motels are an Important tool for testing many hypothosos 
concerning cllmato variability. During recent years a wide variety of 
models have been developed. Complexity of such models ranges between 
the simple energy balance models (e.g. Budyko, 1969; Sellers, 1973) and 
the multi-level primitive equation models (o.g. Manabe ot aK , 1965; 
Kasahara and Washington, 1971; Corby ot aK , 1977; Otto-Blel snor et ajk , 
1982). 

Intermediate complexity models (Klkuchl, 1969; Salmon and 
Hendershott, 1976; Held and Suarez, 1978), with reasonable dynamical end 
physical simplifications, can simulate some aspects of the largest 
scales of atmospheric motion. The computational economy of such models 
provides the opportunity for longer periods of simulation and for more 
extensive testing of physical and dynamical processes. Moreover, such 
models can provide a first insight on atmospheric problems before using 
the complicated general circulation models. Also, Intermediate com- 
plexity models are useful for Interpreting the results of more compli- 
cated models (Chervln, et al_. , 1980). 

In this study a two-level spectral model using pressure as a ver- 
tical coordinate Is developed. The system of equations describing the 
model Is quasl-geostrophlc In linear balance (Lorenz, 1960). The choice 
of global rather than hemispheric model Is due to the fact that the 
latter Is believed to excite anomalous Rossby waves (Roads and 
Somerville, 1982) which could be critical when dealing with climate 
sensitivity studies. 

The physical forcing Is parameterized with reasonable simplicity, to 
include the major forcing mechanisms which develop the large scale 


atmospheric circulation. The solar energy Is specified as a function of 
latitude and time assuming a dally mean zenith angle (Wotherald and 
Manabo, 1972). The amount of solar onergy absorbed by the mod© 1 'fc 
atmosphere and the earth's surface Is calculated using a formula g1«>en 
by Kubota (1972). Longwave radiation forcing of the two layers and the 
surface are calculated using climatological relative humidity and sur- 
face temperature. The mechanical forcifrig of topography Is Introduced in 
the form of a lower boundary vertical velocity. The differential diabe- 
tic heating due to the distribution of land and sea also Is Included. 
The sea surface temperatures are specified using the observed January 
mean values. On continents and Ice surfaces the thermal onergy balance 
equation Is solved for the surface temperature, 
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Fig. 1. Schematic representation of the vertical structure of the model. 
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2. Governing Equations 

The dry flat version of the model structure is basically the same 
as that given by Lorenz (1960) which is a two level, linear balance 
model using pressure as a vertical coordinate. The system of equations 
describing the model retains the nonlinear interactions between depen- 
dent variables. The equations representing the model are the vorticity 
equation, the thermodynamic equation, the thermal wind equation, the 
continuity equation and the water vapor equation. The latter is cal- 
culated at the lower layer only. Static stability is a variable in the 
model's atmosphere and the horizontal wind has both the divergent and 
nondivergent components. 

2.1 Vertical structure of the model (pressure coordinate) 

The model's atmosphere is represented by two levels; 750 mb (£=1) 
and 250 mb (£=3) (Fig. 1). The vertically averaged values are cal- 
culated in the intermediate level 500 mb (£= 2). The lower boundary is 
at the 1000 mb (£=0). 

For a certain level £ the set of equations describing the models 
atmosphere is given by; 


V£ = kx Vi|» £ + V Xr (2.1) 

the vorticity equation 

a 9uj p 

It ^ = - J <v 7 V f > - 7 v 7t * f §r + - (F vV (2 - 2) 

X* x> 

the thermodynamic energy equation 



90 D r 

' J( W" 7 V 7 0jf“;e 8F * ( pf’ 


Q. 

~ + (W. ) +(W ) , (2.3) 

c p n £ v £ 


I 
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the thermal wind equation 

V p o> K ^ <"V- 


(2.4) 


the continuity equation 
0w 

— - + y2 x = o, 

0p x £ 


(2.5) 


and the water vapor equation 

ft ■ ' + E - p c * <Vr 


( 2 . 6 ) 


where = ^ U £’ V 2 ^ the hor ^ zonta1 wind vector, the vertical pres J 

sure velocity, f Is the corlolls parameter, iji Is the stream function, 

Is the velocity potential, Is the potential temperature, q Is the 

water vapor mixing ratio, p £ Is the pressure, pg is the lower boundary 

pressure level (= 1000 mb), P c Is the precipitation rate, E Is the 

surface evaporation rate, Q /c is the diabatic heating rate, c is the 

£ P P 

specific heat at constant pressure, k = R/c p , R is the gas constant, F^, 
W. , S, are the horizontal diffusion of momentum, heat and moisture 
respectively, (F y ) and (W v ) £ are the vertical diffusion of momentum and 
heat, respectively. 

Equations (2.1)-(2.6) are six equations in the 14 unknowns 4^, x £ . 

V V «• V V < F hV < F vV < W hV (W vV <Vr E and P c- The 

evaporation rate, E, is a result of the moisture vertical diffusion from 
the surface while, the precipitation P c is calculated as the excess of 
super saturated moisture In the lower layer. In order to close the set 
(2.1)-(2.6) the diabatic heating and the diffusion terms need to be 
parameterized in terms of the dependent variables. 


2.2 Horizontal diffusion 


From the numerical stability view point the diffusive terms are not 
required when using the spectral method. There Is a requirement to 
Inhibit spurious growth of amplitude at scales close to the point of 
truncation due to spectral blocking (Purl and Bourke, 1974). At a level 
& the horizontal diffusion of momentum, heat and moisture Is param- 
eterized, respectively. 

(f h>» = k k ^ '’’♦j * ! < 2 - 7 > 

<V* = k h 72 V < 2 - 8 > 

CS h ) £ = k h V*q, (2.9) 

where k. Is the lateral eddy diffusion coefficient. The value of k Is 
h n 

+5 -1 

taken to be 1.0x10 m 2 sec (Phillips, 1956). The last term to the 
right side of (2.7) Is due to the effect of spherical earth. 

2.3 Vertical diffusion 

The planetary boundary layer Is a transition layer in the atmo- 
sphere wh'ch separates between the earth surface and the large scale 
atmospheric motions. In this layer, which Is approximately 1 km, the 
fluxes are mainly a consequence of small-scale turbulence and convec- 
tion. In a large scale model It is necessary to utilize the effects of 
the boundary layer to simulate a correct phase and amplitude of the 
ultra-long waves. Parameterized bulk formulas are used here to calcu- 
late the friction dissipation, sensible heat flux and evaporation rate. 
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2.3.1 Parameterization of frictional dissipation 

The two assumptions used for parameterizing the frictional dissipa- 
tion are as follows (Lorenz, 1961): 

a) Surface frictional drag is proportional to the flow In the 
surface layer, 

b) Friction between the two layers Is proportional to the dif- 
ference between the flow of the two layers. 

The friction dissipation, (F )^, Is given by 





( 2 . 10 ) 


where g Is the acceleration of gravity and r Is the rotational stress 

Xr 

at level £.. 

Using the above two assumptions we can have 


r = ^ k V 2 4 j-, 

0 g s ^0 ’ 

and 

r 2 = f 2k d V*«, 3 - *,). 


( 2 . 11 ) 


( 2 . 12 ) 


where Ap(=Pg/2) is the pressure difference between the upper and lower 
levels, and i|» 0 is the surface stream function calcualted by linear 
extrapolation with respect to height (Salmon and Hendershott, 1976). k g 
and 2k d are the coefficients of friction at the underlying surface and 
the surface separating the two layers respectively. k s is given the 

“6 *"1 -7 

value 4x10 sec (Kikuchi, 1969), and k^ is given the value 5x10 
-1 

sec (Charney, 1959). 

Using (2.10), (2.11) and (2.12), and assuming that r^ at the top of 
the atmosphere is equal to zero, we can find the expressions for the 
friction dissipation at the two levels, 
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(F v } = - k s V 2 ip o + 2k d 7 2 (^i 3 - 4^). (2.X3) 

(F y ) = - 2k d V 2 <4i 3 - t| i 1 ), (2-14) 

3 


2.3.2 Parameterization of sensible heat 

Over all surfaces, whether bare land, Ice or water, the vertical 
(turbulent) flux of sensible heat Q g Is determined using the parameteri- 
zation 


Q = p c c . v (T - T ) , 
p d o g a 


(2.15) 


where p is the surface air density, T Is the ground or surface tem- 
s y 

perature (prescribed over the oceans), T fl the surface air temperature, 
c d Is the drag coefficient and jv^jls the absolute value of the surface 
wind. 


The surface air temperature, T , Is extrapolated from the tem- 

a 

perature values at 250 mb and 750 mb with respect to logarithm of the 
pressure level, 


(T a " V/< T a - V " *n(p 0 / Pl )/Jin(p 0 /p 3 ) £ .207 


(2.16) 


The drag coefficient, c rf , is assumed constant taken to be .004 and .001 
over land and water surfaces respectively. By assuming these constant 
values for the drag coefficient we neglected its possible variations 
with the surface wind speed and the terrain height. The absolute value 
of the surface wind, | v q | , Is taken from the rotational part of the 750 
mb wind. A minimum value is specified by 3 m sec to avoid unrealistic 
high surface temperatures (Holloway and Manabe, 1971). 


2.3.3 Parameterization of surface evaporation rate 

The surface evaporation rate, E, Is parameterized In the model as 

6 “ p s c d l u 0 I GW < h »W ' h s q s ( V>’ <2 - 17> 

where q (T ) Is the saturation mixing ratio using the surface tempera- 
s g 

ture, q (T ) the saturation mixing ratio at 1000 mb. The saturation 
s s 

vapor pressure Is calculated using a formula given by Bolton (1980), 
The ground wetness parameter GW is a nondlmensional measure of the 
surface water available for evaporation and varies between 0 and 1. 
Over water and Ice It Is taken as unity, whereas over land surfaces It 
Is taken as .25. The relative humidity In the atmosphere near the 
surface, h , Is given by h = .5 q(T.)/q (T-)+.5, where q(T-) Is the 

mixing ratio in the lower layer. h* Is simply set equal to 1; the 

surface is assumed to be everywhere saturated (the "swamp" lower boun- 
dary condition). 

2.4 Mechanical forcing of topography 

At the top of the model's atmosphere (p=0) the vertical pressure 

velocity tu 4 is taken to be zero. At the lower boundary (1000 mb) io Q 

Introduce t.he mechanical effect of topography, the kinematic condition 

w Q = J(i|/ 1 , P g ), (2.18) 

Is used. Here P g Is the pressure at the terrain height. When computing 

P , the continental elevations smoothed over 5° latitude by 5° longitude 
0 

are used (Berkofsky and Berton, 1955) assuming a standard atmosphere. 
In this relation the advection by the divergent part of the horizontal 
wind is ignored. 

Integration of the continuity equation (2.5) over the depth of the 
model's atmosphere and through its two layers gives the following pres- 


sure velocities 
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ui Q = - ApV a (x x + X 3 ). (2.19) 

^ " f - ? 2 (2 X 3 + Xl ), (2.20) 

and 

u» 3 = - | £ V2 (x 3 ). (2.21) 

It Is convenient to Introduce the new variable x 0 such that 

iu 0 = - Ap V*x 0 . (2.22) 

From (2.19) and (2.22) we get 

X 0 = X 3 + X 3 (2.23) 


The low order truncation used in the model (truncate at either zonal 

wave number 9 or wave number 15) Is considered as a further filter to 

satisfy the quasi-geostrophic approximation, where the vertical velocity 

should be three orders of magnitude less than the horizontal wind > 

(Hal tiner , 1971). 

2. 5 The model 

It is convenient to use as dependent variables the mean potential 
temperature e and the static stability a, the stream functions i|» and t 
for the mean wind and wind shear, so that 0^ ~ 0+cr, 6.^ = 0-o, = i|)+i , 

4^ = i|i~x , x 1 ~ x- Using (2.7) - (2.9), (2.13) and (2.14), the governing 
equations (2.1) - (2.6) become 

~ (V 2 iji) = -J(.|i,V 2 t i.+f)~J(T,V z i)-isV-(fVx 0 )-/v 2 .| )o +k h (V'‘4)+2^), (2.24) 

i 

ii 

i 

A 


^(9 2 X) = -J(>li,7 2 T)-J(t,V 2 i|i+f)+7.(fVx)-HV-(fVx 0 ) + 2 S ^ 0 -2k d V 2 t 


+ k h (V4t + I 2 V2t) * 


(2.25) 


ao 

at 


da 

at 


-J<tt>.())-J<T,(»)+7-<o7x)- J J<Vx 0 *V0+7x 0 'Vo+3oV 2 x 0 ) + k h V 2 0+Q, 
- J(i|* ,a )- J (T , 0 )+Vx • V0-Js(Vx o • V0+Vx o • Va-oV 2 x Q )+k h V 2 o+Q , 


(2.26) 


(2.27) 


§2 

3t = -V‘(kx7(tp-x)+7x)q)+E-P c +k h V 2 q, 


(2.28) 


and 


b c v 2 6 = V.(fVt), 

r 


Ap7 2 x 0 = -J(ijJ“T, P g ), 
i|j 0 = t 


(2.29) 

(2.30) 
(2.3X) 


where 


K K 

b = > 5 [(f) " (|) 3 = . 124 , 


P 0 K P 0 K 

q = Q 3 * C^> Q,]/^ 

is the vertically averaged dlabatic heating per unit mass, and 
P Q K P Q K 

Q 3-^ V /C p 


is the difference in the diabatic heating per unit mass between the two 
layers. 

The above system Is a set of eight equations with eight unknowns 
t, Q, cj, x, X 0 » ^o’ q ‘ This s y stem w '’ 11 be transformed to the spectral 


space using the spherical harmonics as basis functions. 
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3. Thermal Forcing of the Earth-Atmosphere System 

Mechanisms that force the model's atmosphere are either external or 
Internal. The upper layer Is heated by short- and longwave radiation, 
by the lateral diffusion of heat, and by the heat released by a con- 
vective adjustment. The lower layer Is heated by short- and longwave 
radiation, lateral diffusion, sensible heat flux from the surface and by 
latent heat release, and Is cooled by the heat transferred upward by the 
convective adjustment. Evaporation provides a source of water vapor 
which is also diffused and lost through precipitation. 

3,1 Solar radiation 

The Incoming solar radiation at the top of the model's atmosphere 
Is calculated as a function of dally mean zenith angle (Wetherald and 
Manabe, 1972). Diurnal variation of the solar energy is excluded. The 
mean zenith angle z Is given by 

cos z - s 1 ncp sin6 + (co$4> cosfi sin H )/H q , (3.1) 

where $ is the latitude angle, 6 Is the declination angle, and H q Is the 
hour angle given by 

H q = cos ^ ( — tant}» tanfi), (3.2) 

6 = 23.45 sin Zn , (3.3) 

N is the number of days measured from day 0 at 00Z at the first of 
January. 

The Incoming solar radiation at the top of the atmosphere is given 


by 
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S = $ H /n, (3.4) 

oo 0 

S = (f) 2 S c cosz, 4* - fi < \ 

0 4 > - 6 > | 

S Is the solar constant taken to be 1400 w/m z . Recent measurements of 

c 

solar irradlance from earth orbiting satellites (Smith, et al_. , 1982) 

give an average value about 1375 w/m 2 . This value Is about 1.8% less 

than the assumed value. Parameters a and a are the Instantaneous and 

s m 

mean distance of the earth from the sun, respectively, 

-T ~ 1 + -01676 sin 2 n (3.5) 

a 360 

m 

The amount of solar radiation absorbed by the earth's atmosphere system 
Is calculated using a formulae given by Kubota (1972). The solar radi- 
ation absorbed by the atmosphere is given by 

S r - x(l-r a )S 0# , (3.6) 

where x Is the absorptivity of the atmosphere taken to be constant = 
.26. The albedo of the atmosphere, r , is calculated taking Into con- 
slderation the observed mean zonal amount of clouds (Berliand, 1960), 

r = (a + pc)c, (3.7) 

a 

where p Is a constant equal to .38, c is the amount of low and medium 
clouds in tenths of sky cover. Although the model has no explicit 
modulation of the clouds, they are implicitly included through the 
atmospheric albedo which affects the solar energy budget. The parameter 


a is a function of latitude. 


Latitude 


Fig. 2. The zonal average of the calculated solar radiation (w/m 2 ) 
absorbed by the atmosphere (dotted line) and the earth 
(full line) for the first of January- 
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The not solar energy absorbod by the earth's surface to given by 
S s = (1-x) (l-r a ) (l-r s )S M , (3.8) 

whero r Is the January zonal average albedo of the earth's surfaco 

9 

(oceans are not included). The surface albedos are categorized as areas 
of permanent ice (albedo « .8), partial snow in middle and low latitudes 
(albedo = .2 to .3), and dense forests (albedo <= .15). The values of 
different parameters used for the January solar radiative calculation 
are shown in Table 1. 

The above formulae give a global average planetary albedo e 34%. 
Stephens et a_l- , (1981), using satellite observations, estimated the 

global average planetary albedo for January to be 31%. Fig. 2 reveals 
the calculated solar radiation absorbed by the atmosphere and the earths 
surface at the first of January. 


3.2 Longwave radiation 

The calculation of the longwave radiative cooling of the atmosphere 
makes use of a parameterization of the outgoing Infrared radiation 
(Thompson and Warren, 1982). The parameterization comprises clear sky. 
Only two parameters are used to predict cl ear-sky outgoing infrared 
irradiance: surface air temperature (T ) and climatological vertical 

fl 

mean relative humidity (RH). 

The clear sky outgoing infrared irradiance at the top of the 
atmosphere is given by 


L 4 


a +a.T + a„ T 2 + a_T 3 , 
o la 2 a 3 a 


where 


(3.9) 


a „ = b 0n + b l„< RH > + b 2„ (RH > 2 ' 


n = 0, 1, 2, 3. 


(3.10) 


Tho values of the b's are given by, 
b Q0 « 2.34414 x 10 2 , 
b 1Q « -3.47968 x 10 1 , 
b 20 = 1.02790 x 10‘, 
b 01 = 2.60065 x 10°, 
b n = -1.62064 x 10°, 
b 21 = 6.34856 x lo" 1 , 
b Q2 = 4.40272 x 10‘ 3 , 
b 12 = -2.26092 x 10 -2 , 
b 22 = 1.12265 x 10" 2 , 
b Q3 = -2.05237 x 10~ 5 , 
b 13 * -9.670 x 10“ 5 , 
b 23 = 5.62925 x 10" 5 . 

The values of RH used for the January simulation are shown In Table 1. 
These values are interpolated f^nm the values given by Thompson and 
Warren (1982). 

The model's longwave emissivity is divided between the upper and 
lower layer by fraction .4 and .6 respectively. The net longwave ir- 
radiance at the earth's surface (Deardorff, 1978) Is given by 

L 0 - W ■ yB V> < 311 > 

where B is the Stephen Boltzman constant, e is the emissivity of the 

H 

ground surface in the Infrared taken to be equal to .95, and y is the 
parameterization for the effective emissivity of the air which is cal- 


culated from the relation 
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Table 1 


Parameters Used for Solar and Longwave Radiation Calculations 


Parameter 

Latitude 

Clouds 

c 

Atmospheric 
Albedo (r ) 

o 

Surface 

Albedo (r ) 
s 

Average rela- 
tive humidity 
(RH) 

84.1 

0.35 

0.096 

0.8 

.48 

76.5 

0.41 

0.129 

0.8 

. 53 

68.9 

0.48 

0.179 

0.8 

.58 

61.3 

0.54 

0.305 

0.4 

. 6 

53.6 

0.56 

0.343 

0.3 

.59 

45.9 

0.54 

0.316 

0.2 

.58 

38.3 

0.45 

0.248 

0.2 

.54 

30.6 

0.37 

0.185 

0.18 

.46 

23. 

0.28 

0.131 

0.15 

.41 

15.3 

0.29 

0. 145 

0.14 

. 38 

7.7 

0.32 

0.167 

0.14 

.43 

0. 

0.38 

0.207 

0. 14 

.57 

-7.7 

0.36 

0.193 

0.12 

.53 

-15.3 

0.35 

0.183 

0.1 

.46 

-23.0 

0.34 

0.166 

0.1 

.38 

-30.6 

0.36 

0.179 

0.1 

.35 

-38.3 

0.42 

0.227 

0. 1 

.4 

-45.9 

0.51 

0.293 

0.1 

.46 

-53.6 

0.60 

0.377 

0.5 

.50 

-61,3 

0.62 

0.369 

0.5 

.53 

-68.9 

0.55 

0.8 

0.8 

.51 

-76.5 

0.47 

0.8 

0.8 

.46 

-84.1 

0.40 

0.8 

0.8 

.41 
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no 

y = (c + (l-c)x, 67x(1670 q V ), (3.13) 

o 

here the value of c, the cloud fraction, Is assumed as a global average 
equal to .5 and q fl Is the water vapor mixing ratio near the surface. 

3.3 Large scale precipitation and latent heat release 

The model has a moisture content In the lower layer (level 1) only. 
The procedure for large scale precipitation and convective adjustment 
starts after completing each time step of integration. The mixing ratio 
at each grid point of the 750 mb level is examined for super-saturation. 

If q(T.) < yq (T-), then no precipitation or convective adjustment 
takes place. The parameter y represents a specified critical relative 
humidity (y = .85 In this study). T^ Is the temperature at any grid 
point In level 1, and q and q s are the mixing ratio and the saturation 
mixing ratio, respectively. 

On the other hand, If q(T^) > yq s (T^), condensation occurs with the 
associated latent heat release. The temperature T^ will be agumented by 
an increment AT, such that 

AT = (q(T 1 ) - q; (T + AT)), (3.13) 

p 

where q' is the new saturation mixing ratio at the temperature T+AT, 
s 

8q s 

dg = vq s + V AT. (3. 14) 

Using the Cl ausi us-Cl apeyron equation, (3.14) takes the form 

Lq s 

q c = yq s + y AT (3.15) 

where R y is the water vapor gas constant and L Is the latent heat of 
condensation. The rate of condensation (precipitation) per unit mass, 


is given by 
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P c * (q - q;)Mt f < 3 . 16 ) 

where At Is the time step of Integration. Using (3.13), (3.15) and 
(3.16) 


9-vq s L 2 

p ‘*“ <1 + vvi %> 


It Is clear that a relevant form of (3.13) Is 



At. 


(3.17) 


(3.10) 


After the release of latent heat In the lower layer as a result of 

the condensation of water vapor, the atmosphere Is tested to see If 

convective adjustment is required. Convection Is assumed to develop if 

the atmosphere is unstable relative to the moist adiabatic lapse rate 

r , then the temperature of the two levels Is adjusted to stabilize the 
s 

model's atmosphere by cooling the lower layer and warming the upper 
layer, with the vertically averaged temperature conserved. The new 
lapse rate is the same as r . 

3.4 Net heating of the Earth-Atmosphere system 

The way in which the model responds to heating and how It simulates 
the observed atmospheric heat balance are fundamental aspects of Its 
ability to reproduce the seasonal distributions of global climate. From 
the previous discussions we can calculate the different partitions of 
the heating function. 

Of basic importance is the net radiation at the top of the atmo- 
sphere which represents the net gain or loss of both solar and longwave 
radiative energy this may be written as 


N = S - r S 

00 00 a 00 


r s (l - x) (1 - 


r ) S 

a w 


*4‘ 


(3.19) 


On the right side of (3.19) the second and third terms represent the 
amount of solar radiation reflected by the atmosphere and the earth's 
surface, respectively, while the last term is the net outgoing longwave 
radiation at the top of the model atmosphere. 

The net radiation at the earth's surface N s may be Written using 
(3.8) and (3.11) as 

N = (1-x) (1-r ) ( 1-r )S - L . (3.20) 

5 a 5 u 

The net surface heating, B s> Is given by 

8 S = N s - Q s - LE, (3.21) 

It is assumed that B s =0, and the resulting equation Is used to determine 

the surface ground temperature T g . Over the water surfaces, on the 

other hand, the surface temperature is assigned and B is not required 

s 

to be zero. 

The net atmospheric heating may be considered by combining the net 
radiation at the top of the atmosphere (3.19), the net surface heat flux 
(3.21), and the Internal release of latent heat accompanying condensa- 
tion (here precipitation). Recognizing that the surface evaporation 
removes heat from the water source and therefore it is not a part of the 
atmospheric heating, we may write the net heating of the atmosphere, B , 
as 

B = x(l-r ) S + L - L, + Q + LP . (3.22) 

a a °° 0 4 M s c 

This expression for B g is also equal to the sum of the atmospheric 
storage of total energy and the divergence of the atmospheric total 
energy flux. 

Finally, we may combine the net surface heating (3.21) and the net 
atmospheric heating (3.22) in order to get the net heating of the 
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combined earth-atmosphere system. This heating Is given by 

B = N + L(P -E) (3.23) 

ea « c 

This may be regarded as the balance of total energy in the earth- 
atmosphere system. 

3.5 Surface temperature 

The surface temperature, T , is used to calculate the bulk formulae 
(2.15) and (2.17). As mentioned before the surface temperatures of the 
water are specified as the climatological values of January. On land 
and Ice surfaces the temperature is calculated from the surface thermal 
energy balance (3.21) assuming negligible heat capacity of the earth 
(B s =0) (Holloway and Manabe, 1971). Over oceanic locations assumed to 
be covered with ice, B s =0 Is also ussumed, but with a term representing 
the heat conduction through the ice (depending on the difference between 
the Ice surface temperature and the freezing point of water) added to 
the right hand side of (3.21). Over all Ice and snow covered surfaces 
the computed surface temperature Is not permitted to rise above 0°C. In 
such a case the excess heat is assumed to be used In melting. Equation 
(3.21) can take the form 

B a M -Q - LE + I(T - 271.2). (3.24) 

s s s g 

The last term on the right hand side represent the effect of heat con- 
duction from unfrozen water below sea ice in the polar latitudes of the 
Northern Hemisphere. Assuming the thermal conductivity of ice, T = 2.1 

C 

J m 1 °K 1 sec \ the temperature of the underlying water is 271. 2°K and 
the ice layer thickness d = 2 m, then the constant I=T /d=1.05 w/m 2 

C 

°K . This term is needed to prevent unrealistically cold temperatures 
in the Northern Hemisphere polar regions during winter. 
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Appendix I explains the method of solving (3.24). 
iteration method is used and is found to be efficient in 
type of equations. 


The Newton 
solving such 




£ ‘Vi. 


23 


4. Numerical Simulation 

The conventional spectral method Is Galerkin's method based on ex- 
panding the different variables with a truncated series of surface 
spherical harmonics. The method is used for the numerical Integration 
of the hydrodynamlcal equations. Two types of expansion are often used, 
the triangular and rhomboldal truncations. The advantages of the spec- 
tral method over the usual finite difference methods are summarized as 
follows (Machenhauer , 1974): 

1) The nonlinear terms are alias free, which prohibits the exis- 
tence of the nonlinear Instability described by Phillips (1959). 

2) Quadratic area Integral Invariants like the kinetic energy and 
enthalpy also are Invariant for the truncated system, since the 
error fields are orthogonal to the variables. 

3) Linear terms are computed without any truncation error. 

4) No special treatment is required for dealing with the polar 
region when using the vortlcity and divergence fields. By con- 
trast, In the finite difference method the horizontal wind com- 
ponents are discontinuous at the pole. 

5) The friction term of the finite difference methods Is necessary 
to prevent aliasing instability. It also Is necessary for the 
removal of energy from the shortwave end of the spectrum. When 
using the spectral method, it also Is important to prevent blocking 
of energy at the highest wave numbers retained, but in this case 
the purpose is only a simulation of the effect of the small scales 
not retained in the representation. 
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A study by Hoskins and Simmons (1974) compared finite difference 
and spectral models. The study showed that no one method has a super- 
iority in all respects. In comparison with the finite difference model, 
the spectral model gave much improved solutions for the amplitudes and 
phases of the predicted waves. On the other hand, the finite difference 
model gave a more accurate representation of the frontal systems. 

It is of Interest to compare the two types of truncation mentioned 
before, namely the rhomboidal and triangular. For the same zonal wave 
number truncation, the triangular representation has fewer degrees of 
freedom than the rhomboidal and hence requires less computing time. If 
we retain the same degrees of freedom in both the triangular and rhom- 
boidal truncations, the former will be more appropriate for mean zonal 
fields than the latter. At the same time the rhomboidal truncation 
could introduce higher wave numbers, namely the eddies. The same study 
by Hoskins and Simmons (1974) did not give a definite conclusion con- 
cerning the comparison between rhomboidal and triangular truncation. In 
some experiments the rhomboidal truncation gave a more accurate approxi- 
mation to the solution than the triangular truncation. In other experi- 
ments the triangular truncation gave a more efficient description of 
Rossby wave Instability. 

In this study we used the rhomboidal truncation since It gives a 
comparable resolution in both horizontal directions. 

4.1 Spectral method 

The dependent variables »(/, t, x. X Q . 0, o, q are expanded in trun- 
cated series of the form 
M Im I+J 

X(M,X) = \ Z | t 

m=-M n= m j 


(4.1) 




*1 ~ 
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,m 


where K Is any variable being studied, X n are harmonic coefficients, \ 
Is longitude, p Is the sine of latitude, m Is the zonal wavenumber, n Is 
the degree of a spherical harmonic component, n~ Jm| Is a meridional 
wavenumber In the sense that there are n- |m| zero crossings of Y™ 
between equator and pole, M Is the highest zonal wave number retained in 
the truncated series, and J Is the highest value of n- |m | retained In 
the truncated series. Y™ are spherical harmonic functions defined by 


Y m = P m (p) e 1m \ 
n n 


,m 


P n are the Associated Legendre functions of the first kind 

p" (M) . (iSHUi <- ■)' ? ilzy£i" |/Z (M ,.„n 

V M; 1 4n (n+m)! J 2 n n , d(j n+ |m> * 

A spherical harmonic coefficient is defined by 

m * rn 

where Y n is the complex conjugate of Y . 


(4.2) 


(4.3) 


(4.4) 


Y n are orthogonal over the surface of the sphere, l.e. 


Zn +1 


ml 


.... * .... 1 for (m ,n ) = (m.n) 

4tt £ [ 1 Y n V n 1 dfJd * - 0 for (m^n.p / (m,n) * 

and are eigenfunctions of the Laplacian operator 


(4.5) 


. . „ in+ll y". 


(4.6) 


where a is the radius of the sphere. The coefficients for negative and 

positive values of m are related in the following way: 

x -m = m x m* 
n n 


Nonlinear terms are transformed from grid point space to spectral space 
using the full transform method (Machenhauer and Rasmussen, 1972; 
Orszag, 1970). The method Is computationally highly efficient relative 
to the Interaction coefficient method foi* J > 9. 

The procedure for calculating the spectral coefficients of th. 
nonlinear terms using the full transform method Is as follows: 

1) Calculate the nonlinear terms at each grid point In physical 
space. 

2) Transform to the Fourier space at each Gaussian latitude, using 
fast Fourier transform routines, 

3) Transform to the spectral space using the Gaussian quadrature 
formula. 

Highly nonlinear terms, like diabatic heating terms cause problems 
In finding their spectral transforms. This problem Is resolved by using 
the full transform method. They are calculated In physical space, then 
added to the nonlinear dynamic terms, and the whole sum Is transformed 
to spectral space. 

To guarantee an alias-free solution, there are two conditions that 

must be fulfilled (Machenhauer and Rasmussen, 1972). These conditions 

specify the minimum number of zonal grid points, N , and the minimum 

9 

number of Gaussian latitudes, I s> on the sphere: 

N > 3 M + 1 
9 

I s > M + 3/2 J. 

In case of the rhomboidal truncation (M = J) used here, the latter 


condition is 
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For the simulation with wavenumber 9, N « 32 and I = 23. On the other 

u * 

hand, for wave number 15 simulation, N « 40 and I « 40. 

g s 

To transform the system (2.24 - 2.31) to its spectral form, each 
variable is expanded using (4.1). The resulting equations are mul- 

m* 

tiplied by Y n and Integration of both sides is performed using equa- 
tions (4.4 - 4.6). The nonlinear terms are calculated using the trans- 
form method mentioned before. 

The system of equations in its spectral form is given by 


i, m = - 


(-00KV**) - + *n- 


2fim 


, m 


n ' n+Z n m / 

n ( n+l D n+l (x 0 ) n+l 


, n-1 n m, . 

+ — Wn-1> 


(4.7) 


s, .m 

2 l ^0 n 


K h n 


O k 

Lnlll ^ _h m 

a 2 ^n a 2 % 


-*2 


* Hi 3^ < / «■ f , \ ^ / , w-to v % i hi , 2flfn j m 

T n " 7^71) tC-JCc.v^) - J«,.v*t)» n * ^,1 x n 


a / n^*2 r%Hi / vfii , n *** 1 n m / \ 

n l W x 0 , nn * T Wn-l’ 

k 2k 

♦ Aio™ - 2k„ ,'-tn iajli t" . -fi t m 

2 'Mrn d n h a 2 n a 2 n 


(4.8) 


2fi D m x* + D m v m ) 

Vl n+1 x n+l n u n x n-l ; ’ 


0™ = {-J(i|«t0) - J(t,o) - ^(VXQ-Ve + Vx 0 *Vo + 3 oV 2 Xq)}^ 


(4.9) 


" ^ k h 6 n * (v '<° V x»n + ^ • 
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■m 


{- - J(t, 0) - Jj(Vx 0 *VO + VXq'Vct ~ 0 ^ Z Xq) 


m 


+ (Vx-VOT “ n 


Lnill 


(4.10) 


m 




k. cj" + Q" , 
h n M n * 


q m = 

1 n 


(V* ((kxv(iJj-T) t- Vx)q))^ - n • ^ 1 ^ ~ 


k h + < E ' P c>n* 


b c O m = 2Q(^ D m i m + 
p n 'n+1 n+1 T n+1 n 


n m m . 
D n Vl } ' 


(4.12) 


^ o'n 


= *1' 


m 


- - II 

1.6 x 


and 


, .m 

^ x 0^n 


2a 2 

n(n+l) 


u(rt), 


p 

P ’ n ' 
o 


(4.13) 

(4.14) 


where i = V“ 1. The spectral transform of terms of the form V*(f7x) or 
V.(fVx) is shown in Appendix (II). 

It must be noted that by solving (4.0), (4.9) and (4.12) we can 
obtain an equation for x- The equations are simplified and solved as a 
system of tridiagonal matrices (Appendix III) to find the spectral 
coefficients of x that satisfy the linear balance approximation. The 
simplification is needed to treat the term (V*aVx) in (4.9). To do 
this, we split a into its global average [cr], and the deviation from 
this average o' , 

a = [oj + a’ . 

Then 

V*(oVx) = [ct]V 2 x + V*(a'Vx). 


The first term on the right side of the above equation is of a larger 
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order of magnitude and Is added to the other unknown terms, which In- 
clude X’ The smaller, second term, Is considered as a known parameter 
and calculated using the values of x at the previous time step. The 
method is found to be stable. it significantly reduces the number of 
calculations at this stage. 

4.2 Energetics of the model 

The two layer model discussed here conserves the sum of kinetic and 
available potential energy under reversible adiabatic processes (Lorenz, 
1960). If one Introduces the topographical forcing as a lower boundary 
vertical velocity, It Is hard to verify the energy conservation (Burger 
and Rlphagen, 1979). It is only the very simple lower boundary condi- 
tion = 0 (used by Lorenz) at p = 1000 mb that guarantees an energy- 
conserving system. 

The kinetic and available potential energies, KE and AP, respec- 
tively, are expressed in the forms 

KE = ^(Vqj'ViJi + Vt • Vt) (4.15) 

9 


and 


2b c Ap 
AP = E 


[•(O')* ± (a 1 ) 2 ) 




(4.16) 


9 [a] + [ a 2 + (e')2 + (a') 2 ! 

The square brackets [ ] indicate the global area average and the dashes 
indicate the deviation from that average. 

In spectral space the kinetic and available potential energy within 
a spherical harmonic mode are given by 


(KE)?! = (Oi™) 2 + (t™) 2 ) n(n+l) <2-6^), for m > 0 


n ga 


(4.17) 


and 



. 2b c d «°:> 2 +<° n > 2 > ( Mo .) 

0 »" * (a<lll)‘ * (»:>« - (OS)*!’*’ 


for n 1 0 m> 0, where fi nQ = 1 and 6 fl = 0 for m > 0. 


4.3 Initial conditions and time integrations 

The model integration starts from a hypothetical, horizontally 
isothermal, atmosphere at rest with a moist adiabatic lapse rate. The 
model runs for 120 days assuming perpetual solar forcing (first of 
January). This initialization procedure is used in order to reach a 
statistically steady state. After that the solar declination is changed 
daily to simulate the climates of January (days 121-150), February (days 
151-180), and March (days 181-210). These runs are considered as con- 
trol runs for the comparable periods within the experiments. 

The time difference method used is the centered (leap-frog) scheme. 
To avoid the growth of unnecessary computational modes, a time smoother 
was used on the prognostic variables (Asselin, 1972) at every time step. 
The diffusion are calculated using values at the previous time step to 
ensure computational stability. The time step used is 2 hours. Appen- 
dix IV shows a flow diagram of the calculation procedure, 
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5. Summary 

In this report a two-level global spectral model Is' developed. In 
spite of the dynamical and physical simplifications, the model could be 
used to simulate the atmospheric large scale circulation. The model Is 
suitable for climate sensitivity experiments In middle and high lati- 
tudes of both hemispheres. The efflcleht computer runs of the model (30 
day Integration, for wave number 9 truncation, requires about 50 sec of 
CPU time using CRAY-1 machine) enable us to perform many experiments and 
test several hypotheses before using the complicated multilevel primi- 
tive equation models. 

The two levels representing the model's atmosphere are 750 mb and 
250 mb. The surface Is assumed at 1000 mb. The model retains the 
nonlinear Interactions between dependent variables. Nonlinear inter- 
actions are important components of midlatitude synoptic motions. 
Additionally, for climate sensitivity studies nonlinear interactions are 
potentially significant since linear solutions are resonant or nearly 
resonant whiie nonlinear solutions are not. The present model uses a 
moisture budget equation at the 750 mb level with moist convective 
adjustment between the two layers. The advectlon by the divergent wind 
Is retained. Temperature and heat fluxes in each layer can differ 
through a variable static stability. 

The physical forcing Is parameterized with reasonable simplicity to 
Include the major forcing mechanisms which develop the large scale 
atmospheric circulation. The solar energy is specified as a function of 
latitude and time assuming a dally mean zenith angle. Longwave radi- 
ation forcing of the two layers and the surface are calculated. The 
mechanical effects of orography are introduced in the form of a lower 



32 




boundary vertical velocity. The differential dlabatlc heating due to 
the distribution of land and sea also Is Included. The sea surface 
temperatures are specified using the observed January mean values, On 
continents and ice surfaces the thermal energy balance equation Is 
solved for the surface temperature. Both orography and differential 
heating between land and sea are importnat for producing a correct phase 
and amplitude of the middle latitudes ultralong waves In linear atmo- 
spheric models. 

A relatively straightforward extension, not yet attempted, Is the 
parameterization of upper level clouds and their associated radiative 
effects. Such future work is envisaged for studying the role of high 
clouds for short-term climate and the earth's radiation budget. 
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APPENDIX I 

The Solution of the Surface Thermal 
Energy Balance Equation 


Using equations (3.8), (3.11), (3.20), (3.21) and (3.24), the steady- 
state surface thermal energy balance Is represented by 

B s - S s - r g B V * V B V • 0 5 - LE * I (T g - 271.2). 


where 


Q s * Ps C p c dl v oFg - V* 
and LE = L p s c rf |v Q | GW(q g (T ) - h q s (T g )). 


We define 1^ and I such that 


and 


I- - p c c . v_ 
1 s p d 0 


r 2 = L p s c al v ol 


The above equation can be written in the form 


. F( V = S s - 6 g BT g * V B V • I l <T g ' V ' W T g> 

- hq (T ) ) - I(T - 271.2). (Al.l) 

b a g 

(Al.l) Is solved for T , using Newtons iteration method. 

9 

Differentiating (Al.l) with respect to T we obtain 

9 

F ' C V ■ “VV* ■ ■ y^y - < Ai - 2 > 

To calculate the saturation mixing ratio, q (T ) , and its derivative, 

s g 

q'(T ), we use a formula for the saturation vapor pressure, e (Bolton, 
9 s 

1979). This formula provides an accuracy of 0.1% in the range -30°C < 

T < 35°C. 

9 


e (T ) = 6.112 exp (17.67 (T - 273.15)/(T - 29.65)) 

s 9 ,9 g 


(A1.3) 
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.622 0 (T ) 

q (T ) = — 3- 

5 »> P - P S (T 0 ) 


(A1.4) 


Di fferentlatlng (A1.4) with respect to T , 

0 


% pe ; (T q ) 

o' (j \ = — — 2 — SL 

Vg ; CP - e ) 

a 

and using (A1.3), one obtains 


(Al. 5) 


i rT > - 17.67 x 243.15 

s U g j " (T - 29.65) 2 
y 


(Al. 6) 


Substituting (A1.5) into (A1.2), we arrive at 


q.P x 4302.645 

r ' (Tg) = “ 4e g BT g 3 ‘ ! 1 " (p~e (T ))(T~-29. 65) 2 " J * 

** y y 


(Al. 7) 


Using (Al.l) and (A1.7), the solution Is convergent in the form 


v+1 


= T 


F(T V ) 
g 

F* (T V 

g 


(Al. 8) 


where the superscripts v and v+1 indicate successive iteration steps. 
Iteration is performed until F(T ) is less than a small, predetermined 

g 


value. 


APPENDIX II 


Spectral Transform of (V*fVx) 

The term (V*fVx) can be expanded In the form 
V*fVx => Vf*Vx + fv 2 x* 

Since f = 2flp, 

V-fVx = p? (1-M 2 ) §£ + 2npV 2 X . 

If we expand x In terms of spherical harmonics defined by (4.1), 


20 «... m 


0Y 


m 

n 20 


m ..m 


= 75 M x: Cl-M 3 ) rrr ' 75 ll n(n+D x*: mV" 


a* mn n 


3p a z mn 


'n ' n' 


or 


V-(fVx) = 9 X ; {“(n z -l)D^ Y™^ -n(n+2) D’^, VI,). 


a- 1 mn n 

where we have used the two recurrence relations 


»m yfll 

n+1 n+l J 


9Y 


m 


(1 - p., - („♦!> oj - n Y^J. 


,m u m 


and 


wl th 


ap 


p Y m = D m . Y m + D m Y m . , 
M n n+1 n+1 n n-1 


D m = 

n Mn a -r 

Applying the transform operator (4.4) on (A2.3) and using (4.5), 
we obtain 

r- .in ~ 20 . . ^ _ m m f o . . r*m m . 

(V-fVx) n = pp (n(n+2) 0 |>+1 X W1 + (n 2 -l) D„ X„. x ). 


(A2.1) 

then 

(A2.2) 

(A2.3) 

(A2.4) 

(A2.5) 


(A2.6) 


APPENDIX III 


Calculation of the Velocity Potential 

To establish the linear balance approximation, equations (4.8), 
(4.9) and (4.12) need to be solved in order to calculate the array x 
that satisfy the linear balance relation. This appendix describes the 
calculation procedure to find x- 

Using the recurrence formula described In Appendix II, equation 
(4.8) can be written In the form 


, ' \ m 

<o„ * - 


^ / n \ fii /ia / n+2 ^ id m n-1 m > ^ « a « \ 

n(n+l) <R t>n * 20 °n + l x n*l * V D n Vl>' (A3 - 1J 


where (R^)™ the spherical harmonics of the linear and nonlinear terms 
that does not contain x* 

Slmllarily, equation (4.9) can be written in the form 


(0)™ = ( R e )J| + (V-(oVx))^, (A3. 2) 

where (R 0 )^ 1 s the same as the definition of (R^)™ hut for the thermodynamic 
equation. 

The generalized thermal wind equation (4.12) can be differentiated 
with respect to time to give the form 


( 0 ) 


m 

n 


20 .n+2 m •m . n-1 _m • m » 

be ( n+l D n+1 Vl + IT D n t n-l ) ' 
P 


(A3. 3) 


substituting the appropriate indices of (A3.1) and (A3. 2) into (A3. 3), 
we can get the diagnostic equation for x in the form, 


A(n,m)x , ” +2 + B(n,m)x™ 
+ G(n ‘ m)(( Vn-l ! 


+ C(n,m)x|IJ_ 2 + E(n,m)(R i .)™ +1 
(V-(oVx))^ + (R e )™, 


(A3- 4) 
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whom 


D±3 n m n m 

A(ri,m) " be n+1 D n+1 °n+2’ 
P 


■<«.-) ■ 8 1 («S8 <■>*,>* * (n '^ ntl> <£>■>. 


be ' (n+1) 2 vw n+l' 
P 


r , mN _ 4££ rr2 m n m 
C(n,m) - — D n D , 

P 


rt \ 20 n+2 m 

E (n,m) - b(; n+1 D n+1 , 

P 


G(n,m) 


20 n-1 Q m 

be n n 
P 


The system (A3. 4) needs the transformation of (V*oVx) In order to be 
solved. In such case the gaussian elimination method can be used to 
solve for x> However, by making the approximation described in the text 
the system ends to a tridiagonal matrix which is more efficient to solve 
than using the gaussian elimination method. 


APPENDIX IV 
Program Description 


Calculations for this model are contained in three programs. Two 
of them produce input to the model: orography harmonics, ocean tem- 
peratures, legendre polynomial coefficients, gaussian latitudes, gaus- 
sian coefficients. The results of those two programs are stored on the 
files: 


orography harmonics legendre polynomials, etc. 

Wave number 9 ADELH1 ADELH4 

Wave number 15 ADELH2 ADELH3 

The third program calculates the time evolution of the general cir- 
culation. The results of the first 120 days of integration with fixed 
solar radiation for wave number 9 with topography, are stored on file 
ADRES2. The same but without topography is on file ADRES3. Subroutines 
for this program are compiled and stored on file ADELH9 for wave number 
9 and on ADEL15 for wave number 15. 


Time loos. 
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ORIC'r • 

OF POOR QUm.l»v 


Table Al. Schematic representation of the sequence of operations 


heading Input datat geusslan coefficients, sin of gausslan latitudes, 
legendre polynomials and their derivative, ocean temperatures, spectral 
coefficients of topography and results from previous runs (spectral 
coefficients) 

~ """" {■ 


Initial conditions for statistical calculations 

T 


Dally solar forcing * solar radiation absorbed by the earth, 
solar radiation absorbed by atmosphere 


[ 




Mountalr effect * vertical voloclty of the lower boundary (surface 
veloc ity potential) (In spectral space) 


Humidity flux (phys. space) - transform ation to spectral space 

* *" 


I Transformatio n to Fourier s paco _ 


' ]’ 
I 


Transformation to grid points (using FFT) 


Dlabatlc heating terms at gausslan^ latitudes 


| Nonlinear terms at gausslan latitudes 

| Transformation to the Fourier, then to spherical harmonic space 

Adding linear contributions for tendencies (RHS of prognostic 
equations but without terms containing velocity potential) 

■l' 


Solvo for velocity potential to satisfy linear balance 


I Centered time Integration with smoothing" 

T 


Transformation of the vapor mixing ratio and temperature 
to grid point domain 

“ * 


Convective adju s tment 


D 


Transformation of the mixing ratio and temperatures to the 
spectral space _• 

Calculation of x from linear balanco equation, with x ( 1 ,H) 
calculated In sub routine TIME as a bpiiDdnry ^orjd 1 1. 1 op 

i 


|_J. Statistical calculation - monthly ,^zqn nj averaging 


Subprograms 

RDTAPE 


SOLA 

OROG 

ADTO 

TRl 

RM1.RM2 

FFT991 

SURFT.EFAP .SHLT.FLOH 
FFT991 .GUASS 


-► 

S1TER.ADT0 

-K 

TIME 


TR2 

> 

CONVEC 

> 

TRl 

y 

BAL 

-y 

ZAV 


Writing results 



Table A2: The most Important variables In the program. 


Physical Space: 


F i ■ 

f 2 = Is (V^) 
f = a 

r 3 34 

f 4 = Ip ’ 2 * 

F 5 = -d-M 2 ) Sjl 

F 6 " k ^ 

r _ Si 

Y l “ aX 

F 8 * ‘dV> Ip (V a O 
F 9 ■ -d'M 2 ) | 

r - 90 

r 10 " Q\ 

hi = "d"M 2 ) | 

P _ 9 o 
r 12 " 3 K 


F 17 = d-M 2 ) |p (qv) 

F 18 = If “I 1 '’ 

F 19 = 0 
F 20 = 72 « 


F n = q 
'22 


F ?> = -dvi j§ 


p = §11 

h 23 3\ 


F 24 = 


F 25 = 


F = — fi 
r 26 3A 


0 Y 

-<^ z ) 5^ 

8X, 



= -(1-M 2 ) 

= §X 


d\ 


§x 

3p 


where: iJj = stream function 

t = shear 

0 = potential temperature 
a = static stability 
X = velocity potential 
q = mixing ratio 

r) = p /AP = normalized surface pressure 
X = surface velocity potential 
M°= sin£ (£ = latitude) 

K = longitude 


45 






Table A3: The most important variables in the program. 

Spectral Space: 

X — t[j ~ stream function 
TO = t = shear 

PT = o = mean potential temperature 
SI = o = static stability 
Q = q = wake vapor mixing ratio 
RK = x “ velocity potential at 750 mb 
Z = 7 2 ij> 

ZTO = V 2 t 
ZRK = V 2 x 
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Table A4: Catalog of subroutines. 

Subroutines for transformations 


In the following subroutines: 

MM indicates number of points in spectral space in the longitudinal direction 

NN indicates number of points in spectral space in the latitudinal direction 

NG indicates number of points in physical space in the longitudinal direction 

NK indicates number of points in physical space in the latitudinal direction 

SUBROUTINE TR1 (XI ,X,MM,NN,N6,NK) 

TR1 transfers variables from physical to spectral space 
Input: XI (NG , NK) = values in physical space 

Output: X(MM,NN) = spectral coefficients 

SUBROUTINE TR2(X,XI ,MM,NN,N6,NK) 

TR2 transfers variables from spectral to physical space 
Input: X(MM,NN) = spectral coefficients 

Output: X I ( NG , NK ) = values in physical space 

SUBROUTINE GUASS (FMK,FMN,NK,MM,NN) 

GUASS transforms variables of the latitude circles from the Fourier to 
the spherical harmonic domain 

Input: FMK(MM,NK) = Fourier coefficients 

Output: FMN(MM,NN) = spherical harmonics coefficients 

SUBROUTINE RM1(X,K,MM,NN,XM) 

RM1 for given latitude finds ;• ouri er coefficient X(MM) for variable in 
physical space 

Input: X(MM,NN) = variable in spherical harmonic domain 

K = index of latitude 
Output: XM(MM) = Fourier coefficients 

SUBROUTINE RM2(X,K,NN,X,MM) 

RM2 finds Fourier coefficients of the meridional derivative for variable 
X on given latitude 

Input: X(MM,NN) variable in spherical harmonic domain 

Output: XMM(MM) = Fourier coefficient of meridional derivative of X 

SUBROUTINE FFT991(A,W0RK,TRIGS, IFAX, INC, JUMP, N,M, ISIPN) 

FFT991 performs a number of simultaneous real /half-complex Fourier transforms, 
or corresponding inverse transforms. See catalog of NCAR subroutines 
(CRAYLIB library). 
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Subroutines for physical processes 
SUBROUTINE SOLA (NK.ND) 


SOLA calculates solar radiation absorbed by the earth and atmosphere. 

Input: NK = number of gaussian latitudes 

ND = day of year 

Output: QSE = solar radiation absorbed by the earth , . onuMnM/cDcuf'/ 
QSR = solar radiation absorbed by atmosphere) 1 " “MMOM/SRENG/ 


SUBROUTINE 0R0G( X , TO , ETA , XO , MM , NG , NK , ALPH ) 


OROG calculates velocity potential at the surface. 

Input: X(MM,NN) = stream function harmonics 

T0(MM,NN) = shear harmonics 

ETA(MM,NN) = surface pressure divided by pressure increment 
harmonics 

ALPH = parameters regulating height of topography 
Output: X0(MM,NN) = surface velocity potential harmonics 

F24(NG,NK) = laplacian of surface velocity potential 
F25(NG,NK) = meridional derivative of surface velocity potential 
F26(NG,NK) = zonal derivative of surface velocity potential 

SUBROUTINE ADTO (RK,COR,MM,NN) 

ADTO calculates coriolis term with velocity potential 
Input: RK(MM,NN) = velocity potential 

Output: CQR(MM,NN) = V(fVx) 


FUNCTION EVAP(QS,QL1,V1,DRAG) 

EVAP calculates evaporation from surface to the lower layer 

Input.: QS = saturation mixing ratio for surface temperature 

QL1 = saturation mixing ratio in the lower layer of atmosphere 
VI = wind speed in the lower layer of atmosphere 
DRAG = drag coefficient 


SUBROUTINE SURFT(PT1 ,Q1 , VI , K , PTS ,QS , CD , CW, EMS , SFE) 

SURFT calculates surface temperature and saturation mixing ratio for this 
temperature 

Input: PT1 = air temperature at 1000 mb 

Q1 = relative humidity in the lower layer x saturated mixing 
ratio for PT1 

VI = wind speed in the lower layer 
K = latitude index 

PTS = surface temperature from previous time step 

CD = drag coefficient 

CW = wetness parameter 

EMS = surface emissivity of the earth 

SFE = parameter used in calculations of longwave emissivity 
depending on cloud fraction and mixing ratio near the 
surface 

QSE = solar radiation absorbed by the earth 
QSR = solar radiation absorbed by atmosphere 
Output: PTS = surface temperature 

QS = mixing ratio for temperature PTS 
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FUNCTION SNLT ( PTS , PTLS t VI , DRAG) 

SNLT calculates sensible heat flux from the ground to the lower layer 
of atmosphere 

Input: PTS = surface temperature 

PTLS = temperature of the air at 1000 mb 
VI - wind speed in the lower layer 
DRAG = drag coefficient 

FUNCTION FLON (PTS,K) 

FLON calculates clear sky outgoing radiation at the top of the atmosphere 
Input: PTS = surface air temperature 

RH = vertical mean relative humidity (in COMMON/RHLM/ 

K = index of latitude 

SUBROUTINE TIME (X , TO , PT , SI ,Q ,MM,DT ,NTIME) 

TIME makes time step with smoothing 

Input: RHS of eq. 4.7-4.11 (in COMMON/RHS/) 

values of variables from N-l time step (in COMMON/TIMES/) 

DT = time step 
NTIME = number of time step 
MM = max wave number +1 
Output: Values of variables on N+l time step 

X = stream function 
TO =■ shear 

PT = potential temperature 
SI = static stability 
Q = water vapor mixing ratio 

SUBROUTINE BAL(PT,T0,MM,NN,t ) 

BAL calculates shear x from linear balance equation 
Input: PT(MM,NN) = potential temperature 

T0(1,NN) = shear calculated in subroutine TIME 
Output: TQ(MM,NN) = shear satisfying linear balance 

SUBROUTINE SITER(RK t ?.RK,RTO,RPT,GSl,SI ,MM,NG,NK) 

SITER solves equation for velocity potential x in spherical harmonic domain 
Input: RK(MM,NN) = velocity potential from previous time step 

ZRK(MM.NN) = leplacian of velocity potential 
RT0(MM,NN) = R.H.S. of equation for t but without terms 
containing velocity potential 
RPT(MM,NN) = R.H.S. of equation for 0 but without terms 
containing x 

SI(MM,NN) - static stability 
Output: RK(MM.NN) = new value of velocity potential 
ZRK(MM,NN) = new value of V 2 x 

G$1(MM,NN) = velocity potential term in equation for 0 
F13(NG,NK) = meridional derivative of x 
F14(NG,NK) = zonal derivative of x 
F20(NG,NK) = V 2 x 



SUBROUTINE ZAV( RTT , NG ,NK , AV) 


ZAV calculates zonal average of variable RTT 

Input: NTT ( NG , KK) £ variable In physical space 

Output: AV(NK) = zonal average of RTT 

SUBROUTINE CONVEC(QG , PTG , SIG , PROP , NG , NK.TIM) 

CONVEC makes convective adjustment and calculates precipitation rate 
Input: QG(NG,NK) = mixing ratio before convective adjustment 

PTG(NG,NK) = potential temperature at 500 mb before convective 
adjustment 

SIG(NG.NK) = static stability before convective adjustment 
TIM » time step 

Output; PTG(NG.NK) = potential temperature at 500 mb after convective 

adjustment 

SIG(NG,NK) = static stability after convective adjustment 
PRCP(NG,NK) - precipitation rate 


